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We use determinant Quantum Monte Carlo simulations and exact diagonalization to explore 
insulating behavior in the Hubbard model with a bimodal distribution of randomly positioned local 
site energies. From the temperature dependence of the compressibility and conductivity, we show 
that gapped, incompressible Mott insulating phases exist away from half filling when the variance 
of the local site energies is sufficiently large. The compressible regions around this Mott phase are 
metallic only if the density of sites with the corresponding energy exceeds the percolation threshold, 
but are Anderson insulators otherwise. 
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Introduction 

The translationally invariant Hubbard model has long 
been studied as a model of itinerant magnetism and 
(Mott) insulating behavior. More recently, the possibility 
of unconventional [d— wave) superconductivity and spon- 
taneously occurring charge inhomogeneities (stripes and 
checkerboards) has been explored, especially in the con- 
text of high temperature superconductivity^ Including 
disorder in the Hubbard Hamiltonian, for example in the 
form of a distribution of bond or site energies, imposes 
charge inhomogeneity externally and allows for the explo- 
ration of a number of other interesting phenomena such 
as the formation of Anderson insulating phases, possible 
transitions to metallic behavior driven by interactions^, 
and the influence of disorder on magnetic correlations 4 :. 

A particularly interesting suggestion made recently^ 
concerns the possibility of Mott insulating phases away 
from half-filling in a Hubbard model corresponding to a 
binary alloy, that is, for which the probability distribu- 
tion of site energies is bimodal, P(e^) = x5(ti + A/2) + 
(1 — x)5(ei — A/2). The idea is that if A is sufficiently 
large compared to the bandwidth W, the non-interacting 
density of states will be split and an insulating gap will 
separate two density of states peaks of weight x and 
1 — x. When an on-site repulsion U is turned on, these 
two peaks may in turn be Hubbard split by U. Thus 
in the limit A > U > W one can have Mott insulating 
phases at incommensurate densities p — x and p = 1 + x 
which correspond to half-filling the two alloy subbands. 
A further interesting aspect of these Mott insulators is 
that they likely occur in the absence of antiferromag- 
netic ordering and its associated symmetry breaking, a 
phenomenon which complicates the metal-insulator tran- 
sition in the translationally invariant Hubbard model at 
p = 1. Possible experimental realizations of the binary al- 
loy Hubbard Hamiltonian in two dimensions include Co- 
Fe monolayers^. The nature of magnetism in such sys- 
tems has been explored by first principles calculations 8 . 

The previous studies of Mott transitions off half-filling 
within tight binding models were with dynamical mean 
field theory (DMFT). In this letter we will re-examine 
the physics of the binary alloy Hubbard model using de- 



terminant Quantum Monte Carlo (DQMC) simulations 
and exact diagonalization. While these methods are re- 
stricted to finite size lattices, they allow us to examine 
some of the aspects of the effects of randomness like An- 
derson localization which are not accessible with DMFT. 
The specific Hamiltonian we study is: 

(ij)<r i 

Here cj a (ci a ) are the usual fermion creation (destruction) 

operators for spin a on site I, ni a — cj^cia- is the num- 
ber operator, and refers to near neighbor pairs on a 
two dimensional square lattice, t, p and U are the elec- 
tron hopping, chemical potential, and on-site interaction 
strength, respectively, and e; is a local site energy given 
by the bimodal distribution described previously. The 
bandwidth is W = 8t when A = U = 0. 

This paper is organized as follows: We will first de- 
scribe some of the details of our computational method- 
ology. We then show results for the density as a func- 
tion of chemical potential which illustrate the appear- 
ance of Mott plateaus off half filling, and also demon- 
strate the consistency of DQMC and direct diagonal- 
ization. Results for the participation ratio in the non- 
interacting limit suggest that the Mott plateaus at p = x 
and p = 1 + x could in fact be rather different, a conclu- 
sion which we then confirm by calculating the tempera- 
ture dependence of the conductivity. Finally, we exam- 
ine the critical hopping t c required to destroy the Mott 
plateau. We conclude by constructing the corresponding 
phase diagram. 



Computational Methods 

We study the alloy Hubbard Hamiltonian with exact 
diagonalization and DQMCS. The former approach is a 
standard application of the Lanczos algorithm to deter- 
mine exactly the ground state wave function. We use 
N = 8 site lattices. In this case N is sufficiently small 
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that we can sum over all disorder realizations. In order 
to reduce finite size effects, we employ the boundary con- 
dition averaging method^. In the noninteracting limit, 
averaging over different hopping phases at the bound- 
ary for a finite lattice reproduces the thermodynamic 
limit spectrum exactly. For U nonzero, the finite size 
effects, while not eliminated, are dramatically reduced. 
Specifically, we implement a 2x4 cluster with two bound- 
ary phases [<j> x and <fi y ), one boundary phase for each of 
the two independent, orthogonal boundaries. An average 
over the phase space area encompassed by <p x = {0, 2tt} 
and <f)y = {0, 2ir} was done by selecting 10 to 100 {4> x , (f> y ) 
pairs. 

In the DQMC approach, the partition function Z is 
expressed as a path integral by discretizing the inverse 
temperature /?. The on-site interaction is then replaced 
by a sum over a discrete Hubbard-Stratonovich field 11 . 
The resulting quadratic form in the fermion operators 
can be integrated out analytically, leaving an expression 
for Z in terms of a sum over all values of the Hubbard- 
Stratonovich field with a summand (Boltzmann weight) 
which is the product of the determinants of two matrices 
(one for spin up and one for spin down). The sum is 
sampled stochastically using the Metropolis algorithm. 
We present results for 6 x 6 lattices. We average over 
5 — 10 realizations of the local site energies. 

Equal time operators such as the density and energy 
are measured by accumulating appropriate elements, and 
products of elements, of the inverse of the matrix whose 
determinant gives the Boltzmann weight. For the con- 
ductivity, <7dc, we employ an approximate procedure^, 
which allows <7d c to be computed from the wavevector q- 
and imaginary time r-dependent current-current correla- 
tion function A xx (q, r) without the necessity of perform- 
ing an analytic continuation^, 

/3 2 

<7dc = — A ra (q = 0,T = /3/2) . 

Here (3 = 1/T, A xx (q,r) = (j x {q, r)j x (-q, 0)), and 
jx (q, T ) the q, r-dependent current in the x-direction, is 
the Fourier transform of, 

This approach has been extensively tested for the 
superconducting-insulator transition in the attractive 
Hubbard model 1 ^, as well as for metal-insulator tran- 
sitions in the repulsive model2*i^. 

In order to get further insight into the physics of An- 
derson localization, we also diagonalize the noninteract- 
ing system on lattices as large as N — 64 x 64. We 
characterize the properties of the noninteracting eigen- 
functions \<p n ) through the scaling of the participation 
ratio, 

N _ x 




-4-3-2-1 1 2 3 4 



FIG. 1: (Color online) Density as a function of chemical po- 
tential for W = 0.8, U = 2, A = 4 and x = 0.25. Solid line: 
ground state exact diagonalization of an eight site cluster; 
and dashed line: DQMC on 6 x 6 clusters with temperature 
T = 1/16. In the former case, results are averaged over all 
configurations with two sites with €a = —2 and six sites with 
€b — +2, and over different choices of the boundary phases. 
In the latter case, results are given for a single realization. 
The two methods give very similar results, with the DQMC 
somewhat rounded by finite temperature. 



For an eigenfunction perfectly localized at a site io, 
( i \ 4> n ) = &i.i a we have, V n = lj while for a perfectly delo- 
calized eigenfunction (i \<p n ) = 1/yjV, we have V n = N. 
In general, V n is a measure of the extent of the eigen- 
function, that is, the number of sites for which ( i \<p n ) is 
non-negligible. 



Results 

We first demonstrate the existence of the Mott and 
band insulating phases by looking at the density as a 
function of chemical potential. Fig. 1 shows the case 
U = 2, A = 4 and x = 0.25. We see a Mott insulating 
plateau extending from // = (—A — U)/2 = —3 to /i = 
(—A + U)/2 = —1, in which the total density is p — 
x = 0.25. At (i = (-A + U)/2 = -1, the lower alloy 
subband becomes doubly occupied and p = 2x = 0.5. A 
second plateau then reflects the 'band' gap which must be 
surpassed to begin occupying the upper alloy subband, 
which is, like the lower subband, also Hubbard split. 

In Fig. 1, the Mott plateaus revealed in p(fi) at p = 
x = 0.25 and p = 1 + x = 1.25 appear rather similar. 
We now argue that the nature of the states is, instead, 
quite different. For a square lattice in d = 2 the per- 
colation threshold is x c = 0.5928 15 . We therefore might 
expect that the non-interacting states with site energy 
ei = —A/2 = —2, out of which the p = x plateau is 
built, are localized, since the constitute only a fraction 
x = 0.25 < x c of the sites in the lattice. Meanwhile, the 
non-interacting states with site energy e; = A/2 = +2, 
out of which the p = 1 + x plateau is built, are delo- 
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FIG. 2: (Color online) Upper panel: Participation ratio P as 
a function of the eigenenergy E for the alloy Hubbard Hamil- 
tonian with U — 0, A = 4, and t = 1. Results are shown for 
lattice sizes varying from N — 16 x 16 to N = 64 x 64. P(E) 
is small for the lower alloy band which has only x = 0.25 < x c 
of the sites, but P(E) in the upper alloy band has a signifi- 
cant fraction of TV. The inset shows that P(E)/N scales to a 
non-zero value as N — > oo for the upper band and zero for the 
lower band. Lower panel: The same for x = 0.45. Here both 
alloy subbands have a density below the percolation thresh- 
old, and both participation ratios scale to zero. 



calized. This is illustrated in Fig. 2 where we show the 
participation ratio of the noninteracting system. States 
with energies corresponding to the upper alloy subband, 
which has a density exceeding the percolation threshold, 
extend over a macroscopic portion of the lattice. States 
in the lower alloy subband are localized. 

We now use the temperature dependence of the con- 
ductivity to argue that the distinction between the two 
Mott insulating plateaus is preserved when the interac- 
tion U is turned on. Fig. 3 gives Udc as a function of fj, 
for three different temperatures T = 1/8, T = 1/12 and 
T = 1/16. The conductivity is zero in both the Mott and 
band insulator phases. In the regions bracketing the up- 
per alloy band, Ode is relatively large, and increases as T 
is lowered. That is, these regions are metallic. For den- 
sities bracketing the lower alloy subband, <7dc is a factor 
of four smaller, and increases much less noticeably as T 
is lowered. Because of the sign problem, we are not able 




FIG. 3: (Color online)The conductivity <7dc is shown as a 
function of chemical potential for several temperatures. Pa- 
rameter values are as in Fig. 1. The upper alloy band, whose 
sites have a density larger than the percolation threshold, has 
a large Od c , which also increases as T is reduced (inset). The 
insulating phases have a^ c ~ 0. The conductivity near the 
lower alloy subband is much smaller, and much more weakly 
temperature dependent than the upper. 
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FIG. 4: (Color online)The length of the insulating plateaus 
Afi is shown as a function of hopping t for U — 2, A = 4 
and x — 0.25, using exact diagonalization of TV = 8 site clus- 
ters using boundary condition averaging. The Mott insulator 
plateaus at p — x and p — 1 + x are more robust than the 
band insulator plateau at p = 2x. 



to obtain data for lower T. However, we believe that, 
as in the case of the Hubbard model with random site 
or bond energie a 2 i 14 , the conductivity will turn over and 
decrease as T is lowered further, reflecting the insulating 
character of the states. 

We have presented DQMC results at U = 2, A = 4 
which correspond to on-site interaction and alloy site en- 
ergy separation about twice and four times the band- 
width, respectively. The reason for these strong coupling 
values is that for larger t it is not possible to reach low 
enough temperatures to see clear plateaus in the den- 
sity versus chemical potential plots. However, we argued 
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FIG. 5: Phase diagram of U/t vs. A/t. I: Band Insulator at 
p — 0.50 and no Mott Insulator. II: Metallic phase. Ill: Band 
Insulator at p — 0.50 and Mott Insulators at p — 0.25 and 
1.25. IV: Mott Insulator at p = 0.25 and 1.25 and no Band 
Insulator. V: U > A, band insulator at p = 0.25 and 1.25 
and Mott insulator at p — 1.00. 

that diagonalization and DQMC gave consistent results 
(see Fig. 1), and we will now use the former approach 
to generate the ground state phase diagram at weaker 
couplings. As before, we average over different boundary 
condition phases to reduce finite size effects. 

Fig. 4 shows the length of the three plateaus as the 
hopping t is increased for U = 2, A = 4, and x = 0.25. 
The two Mott plateaus appear to vanish at roughly the 
same hopping strength, t w 1. The band insulator is less 
robust, and is destroyed when quantum fluctuations are 
only about half as strong. The reason is that when hop- 
ping off of one of the lower energy alloy sites for a Mott 
insulator you have to pay a cost of A = 4. However, when 
you hop off of one of these sites for the band insulator 
you only have to pay the cost of A — U — 2. This greater 
ease of such charge fluctuations for the band insulator 
makes its destruction by increasing t occur earlier. 

Similar data for other choices of U and A allow us to 
generate the ground state phase diagram. In Fig. 5 we 
show four possible phases for the region of which U < A. 
As in Fig. 4, large t that correspond to small U/t and 
A/t will suppress the insulating plateaus which results 
in a metallic phase. Taking the limiting case for small 
U/t and large A/t, we find one band insulating plateau 



at p = 2x. Inversely, for small A/t and large U/t, we 
find two Mott phases at p = x and p = 1 + x. For 
large U/t and A/t, both Band and Mott insulator phases 
coexist. The case for U > A will correspond to two band 
insulators at x and 1 + x and a Mott insulator at half- 
filling for sufficiently large values of U/t and A/t. 



Conclusions 

In this paper we have presented DQMC and diago- 
nalization results for the phase diagram of the Hubbard 
model with binary alloy disorder. In agreement with 
previous treatments 5 '-, we find Mott insulating behavior 
away from half-filling when the separation of the two site 
energies exceeds U. We extended the earlier results to 
characterize the nature of the compressible states above 
and below the Mott insulating plateaus by showing that 
their conductivity markedly differs. Together with the re- 
sults for the participation ratio, our data suggest that for 
x < x c the lowest Mott plateau separates two compress- 
ible Anderson insulating regions, while the upper Mott 
plateau separates two compressible metallic phases. 

Our results have focused primarily on x — 0.25, but 
different behavior would emerge for other values of a»i&. 
For example, choosing a value of x — 0.70 > x c would 
not only produce Mott plateaus at p = 0.70 and 1.70, 
but also create a lower Mott gap that are surrounded 
by metallic phases. Consequently, the upper Mott gap 
(1 — x = 0.30 < x c ) would be in between two Anderson 
insulating states. 

A closer examination of magnetic correlations in this 
model is of interest, and will be the subject of future 
work. While we expect that disorder will suppress mag- 
netism, as will the fact that the Mott phases are away 
from commensurate fillings, it is also the case that disor- 
der can increase the exchange constant J and hence the 
Neel temperature in certain circumstances^^. Exam- 
ining the real space magnetic correlations using DQMC 
would be a useful complement to previous DMFT stud- 
ies. 
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